import pandas as pd
import geopandas as gpd
grid_path=r"C:\Users\User\Documents\EHI\Projects\VAR Modelling\Data Mgmt\500 grids\grid_500m_clipped.shp"
grid_layer = gpd.read_file(grid_path)
grid_layer.head()
| OBJECTID | Id | Shape_Leng | Shape_Area | geometry | |
|---|---|---|---|---|---|
| 0 | 1 | 0 | 2000.000000 | 250000.000000 | POLYGON ((18162.733 31755.775, 17662.733 31755... |
| 1 | 2 | 0 | 2000.000000 | 250000.000000 | POLYGON ((23162.733 42755.775, 23662.733 42755... |
| 2 | 3 | 0 | 2000.000000 | 250000.000000 | POLYGON ((37662.733 40255.775, 37162.733 40255... |
| 3 | 4 | 0 | 2961.516312 | 226945.081982 | MULTIPOLYGON (((28662.733 48255.775, 28235.903... |
| 4 | 5 | 0 | 1974.913982 | 248059.096956 | POLYGON ((29162.733 28255.775, 29162.733 27755... |
grid_layer=grid_layer.set_crs(epsg="3414")
grid_layer.info()
<class 'geopandas.geodataframe.GeoDataFrame'> RangeIndex: 3499 entries, 0 to 3498 Data columns (total 5 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 OBJECTID 3499 non-null int64 1 Id 3499 non-null int64 2 Shape_Leng 3499 non-null float64 3 Shape_Area 3499 non-null float64 4 geometry 3499 non-null geometry dtypes: float64(2), geometry(1), int64(2) memory usage: 136.8 KB
df = pd.read_csv("DEN-Hexagon DBF 2013_2020_complete.csv")
df.head()
| CASE_ID | CLUSTER_ID | Onset_EYear | Onset_EWeek | Clean_Onset | ONSET_DATE | NEA_ONSET_ | NOTIFY_DAT | POSTAL | X output | Y output | cell_id | friendly_n | Shape_Leng | Shape_Area | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | DENF-2013-000025 | 0 | 2012 | 51 | 17/12/2012 | 17/12/2012 | 17/12/2012 | 31/12/2012 | 650536 | 18857.18765 | 37512.04877 | 1334 | HONG KAH NORTH | 0.016238 | 0.000019 |
| 1 | DENF-2013-000018 | 20130008 | 2012 | 51 | 21/12/2012 | 31/12/2012 | 21/12/2012 | 31/12/2012 | 536467 | 32943.97740 | 36955.30429 | 2826 | UPPER PAYA LEBAR | 0.016238 | 0.000019 |
| 2 | DENF-2013-000048 | 0 | 2012 | 51 | 22/12/2012 | 30/12/2012 | 22/12/2012 | 3/1/2013 | 550421 | 32663.67780 | 36774.03736 | 2826 | UPPER PAYA LEBAR | 0.016238 | 0.000019 |
| 3 | DENF-2013-000109 | 0 | 2012 | 51 | 22/12/2012 | 3/1/2013 | 22/12/2012 | 4/1/2013 | 550422 | 32750.05966 | 36774.55328 | 2826 | UPPER PAYA LEBAR | 0.016238 | 0.000019 |
| 4 | DENF-2013-000017 | 0 | 2012 | 51 | 20/12/2012 | 28/12/2012 | 20/12/2012 | 1/1/2013 | 534969 | 33192.61720 | 36691.75100 | 2861 | TAI SENG | 0.016238 | 0.000019 |
gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df["X output"], df["Y output"]))
gdf=gdf.set_crs(epsg="3414")
gdf.head()
| CASE_ID | CLUSTER_ID | Onset_EYear | Onset_EWeek | Clean_Onset | ONSET_DATE | NEA_ONSET_ | NOTIFY_DAT | POSTAL | X output | Y output | cell_id | friendly_n | Shape_Leng | Shape_Area | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | DENF-2013-000025 | 0 | 2012 | 51 | 17/12/2012 | 17/12/2012 | 17/12/2012 | 31/12/2012 | 650536 | 18857.18765 | 37512.04877 | 1334 | HONG KAH NORTH | 0.016238 | 0.000019 | POINT (18857.188 37512.049) |
| 1 | DENF-2013-000018 | 20130008 | 2012 | 51 | 21/12/2012 | 31/12/2012 | 21/12/2012 | 31/12/2012 | 536467 | 32943.97740 | 36955.30429 | 2826 | UPPER PAYA LEBAR | 0.016238 | 0.000019 | POINT (32943.977 36955.304) |
| 2 | DENF-2013-000048 | 0 | 2012 | 51 | 22/12/2012 | 30/12/2012 | 22/12/2012 | 3/1/2013 | 550421 | 32663.67780 | 36774.03736 | 2826 | UPPER PAYA LEBAR | 0.016238 | 0.000019 | POINT (32663.678 36774.037) |
| 3 | DENF-2013-000109 | 0 | 2012 | 51 | 22/12/2012 | 3/1/2013 | 22/12/2012 | 4/1/2013 | 550422 | 32750.05966 | 36774.55328 | 2826 | UPPER PAYA LEBAR | 0.016238 | 0.000019 | POINT (32750.060 36774.553) |
| 4 | DENF-2013-000017 | 0 | 2012 | 51 | 20/12/2012 | 28/12/2012 | 20/12/2012 | 1/1/2013 | 534969 | 33192.61720 | 36691.75100 | 2861 | TAI SENG | 0.016238 | 0.000019 | POINT (33192.617 36691.751) |
gdf_latlong = gdf.to_crs(epsg = "4757")
gdf_latlong.head()
| CASE_ID | CLUSTER_ID | Onset_EYear | Onset_EWeek | Clean_Onset | ONSET_DATE | NEA_ONSET_ | NOTIFY_DAT | POSTAL | X output | Y output | cell_id | friendly_n | Shape_Leng | Shape_Area | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | DENF-2013-000025 | 0 | 2012 | 51 | 17/12/2012 | 17/12/2012 | 17/12/2012 | 31/12/2012 | 650536 | 18857.18765 | 37512.04877 | 1334 | HONG KAH NORTH | 0.016238 | 0.000019 | POINT (103.75116 1.35552) |
| 1 | DENF-2013-000018 | 20130008 | 2012 | 51 | 21/12/2012 | 31/12/2012 | 21/12/2012 | 31/12/2012 | 536467 | 32943.97740 | 36955.30429 | 2826 | UPPER PAYA LEBAR | 0.016238 | 0.000019 | POINT (103.87774 1.35048) |
| 2 | DENF-2013-000048 | 0 | 2012 | 51 | 22/12/2012 | 30/12/2012 | 22/12/2012 | 3/1/2013 | 550421 | 32663.67780 | 36774.03736 | 2826 | UPPER PAYA LEBAR | 0.016238 | 0.000019 | POINT (103.87522 1.34885) |
| 3 | DENF-2013-000109 | 0 | 2012 | 51 | 22/12/2012 | 3/1/2013 | 22/12/2012 | 4/1/2013 | 550422 | 32750.05966 | 36774.55328 | 2826 | UPPER PAYA LEBAR | 0.016238 | 0.000019 | POINT (103.87600 1.34885) |
| 4 | DENF-2013-000017 | 0 | 2012 | 51 | 20/12/2012 | 28/12/2012 | 20/12/2012 | 1/1/2013 | 534969 | 33192.61720 | 36691.75100 | 2861 | TAI SENG | 0.016238 | 0.000019 | POINT (103.87998 1.34810) |
hexagon_path = r"C:\Users\User\Documents\EHI\Projects\VAR Modelling\Data Mgmt\hexagons\hexagons.shp"
hexagons_layer = gpd.read_file(hexagon_path)
hexagons_layer.crs
<Geographic 2D CRS: EPSG:4757> Name: SVY21 Axis Info [ellipsoidal]: - Lat[north]: Geodetic latitude (degree) - Lon[east]: Geodetic longitude (degree) Area of Use: - name: Singapore - bounds: (103.59, 1.13, 104.07, 1.47) Datum: SVY21 - Ellipsoid: WGS 84 - Prime Meridian: Greenwich
hexagons_layer.info()
<class 'geopandas.geodataframe.GeoDataFrame'> RangeIndex: 4083 entries, 0 to 4082 Data columns (total 11 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 OBJECTID 4083 non-null int64 1 TARGET_FID 4083 non-null int64 2 centroid 4083 non-null object 3 cell_id 4083 non-null object 4 friendly_n 4083 non-null object 5 Shape_Leng 4083 non-null float64 6 Shape_Area 4083 non-null float64 7 cell_id_1 4083 non-null int64 8 X 4082 non-null float64 9 Y 4082 non-null float64 10 geometry 4082 non-null geometry dtypes: float64(4), geometry(1), int64(3), object(3) memory usage: 351.0+ KB
hexagons_layer.head()
| OBJECTID | TARGET_FID | centroid | cell_id | friendly_n | Shape_Leng | Shape_Area | cell_id_1 | X | Y | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 0 | { "type": "Point", "coordinates": [ 103.674855... | 463 | GUL BASIN | 0.016238 | 0.000019 | 463 | 103.675 | 0.0 | POLYGON ((103.67216 1.29272, 103.67351 1.29507... |
| 1 | 2 | 1 | { "type": "Point", "coordinates": [ 103.678899... | 514 | SAFTI | 0.016238 | 0.000019 | 514 | 103.679 | 0.0 | POLYGON ((103.67620 1.33270, 103.67755 1.33505... |
| 2 | 3 | 2 | { "type": "Point", "coordinates": [ 103.691040... | 634 | SHIPYARD | 0.016238 | 0.000019 | 634 | 103.691 | 0.0 | POLYGON ((103.68834 1.30213, 103.68969 1.30448... |
| 3 | 4 | 3 | { "type": "Point", "coordinates": [ 103.840759... | 2410 | OXLEY | 0.016238 | 0.000019 | 2410 | 103.841 | 0.0 | POLYGON ((103.83806 1.29508, 103.83941 1.29743... |
| 4 | 5 | 4 | { "type": "Point", "coordinates": [ 104.018803... | 3717 | CHANGI BAY | 0.016238 | 0.000019 | 3717 | 104.019 | 0.0 | POLYGON ((104.01611 1.31388, 104.01745 1.31623... |
hexagons_layer.dropna(inplace=True)
hexagons_layer["cell_id"]=hexagons_layer["cell_id"].astype(int)
hexagons_layer=hexagons_layer[hexagons_layer["cell_id"].isin(df["cell_id"])]
hexagons_layer.info()
<class 'geopandas.geodataframe.GeoDataFrame'> Int64Index: 1563 entries, 3 to 3742 Data columns (total 11 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 OBJECTID 1563 non-null int64 1 TARGET_FID 1563 non-null int64 2 centroid 1563 non-null object 3 cell_id 1563 non-null int32 4 friendly_n 1563 non-null object 5 Shape_Leng 1563 non-null float64 6 Shape_Area 1563 non-null float64 7 cell_id_1 1563 non-null int64 8 X 1563 non-null float64 9 Y 1563 non-null float64 10 geometry 1563 non-null geometry dtypes: float64(4), geometry(1), int32(1), int64(3), object(2) memory usage: 140.4+ KB
hexagons_layer.to_crs(epsg="3414")["geometry"].length
3 1801.332902
20 1801.332975
32 1801.332740
34 1801.332770
76 1801.332915
...
3626 1801.332947
3644 1801.332834
3648 1801.332858
3669 1801.332866
3742 1801.332902
Length: 1563, dtype: float64
hexagons_layer.to_crs(epsg="3414")["geometry"].area/ 10**6
3 0.234173
20 0.234173
32 0.234173
34 0.234173
76 0.234173
...
3626 0.234173
3644 0.234173
3648 0.234173
3669 0.234173
3742 0.234173
Length: 1563, dtype: float64
import folium
from folium.plugins import HeatMap
map_folium = folium.Map([1.35255,103.82580],zoom_start=11.4,control_scale=True)
for _, r in hexagons_layer.iterrows():
#without simplifying the representation of each borough, the map might not be displayed
#sim_geo = gpd.GeoSeries(r['geometry'])
sim_geo = gpd.GeoSeries(r['geometry']).simplify(tolerance=0.001)
geo_j = sim_geo.to_json()
geo_j = folium.GeoJson(data=geo_j,style_function=lambda x: {'fillColor': 'orange'})
# folium.Popup(r['BoroName']).add_to(geo_j)
geo_j.add_to(map_folium)
map_folium
gdf.head()
| CASE_ID | CLUSTER_ID | Onset_EYear | Onset_EWeek | Clean_Onset | ONSET_DATE | NEA_ONSET_ | NOTIFY_DAT | POSTAL | X output | Y output | cell_id | friendly_n | Shape_Leng | Shape_Area | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | DENF-2013-000025 | 0 | 2012 | 51 | 17/12/2012 | 17/12/2012 | 17/12/2012 | 31/12/2012 | 650536 | 18857.18765 | 37512.04877 | 1334 | HONG KAH NORTH | 0.016238 | 0.000019 | POINT (18857.188 37512.049) |
| 1 | DENF-2013-000018 | 20130008 | 2012 | 51 | 21/12/2012 | 31/12/2012 | 21/12/2012 | 31/12/2012 | 536467 | 32943.97740 | 36955.30429 | 2826 | UPPER PAYA LEBAR | 0.016238 | 0.000019 | POINT (32943.977 36955.304) |
| 2 | DENF-2013-000048 | 0 | 2012 | 51 | 22/12/2012 | 30/12/2012 | 22/12/2012 | 3/1/2013 | 550421 | 32663.67780 | 36774.03736 | 2826 | UPPER PAYA LEBAR | 0.016238 | 0.000019 | POINT (32663.678 36774.037) |
| 3 | DENF-2013-000109 | 0 | 2012 | 51 | 22/12/2012 | 3/1/2013 | 22/12/2012 | 4/1/2013 | 550422 | 32750.05966 | 36774.55328 | 2826 | UPPER PAYA LEBAR | 0.016238 | 0.000019 | POINT (32750.060 36774.553) |
| 4 | DENF-2013-000017 | 0 | 2012 | 51 | 20/12/2012 | 28/12/2012 | 20/12/2012 | 1/1/2013 | 534969 | 33192.61720 | 36691.75100 | 2861 | TAI SENG | 0.016238 | 0.000019 | POINT (33192.617 36691.751) |
grid_layer.head()
| OBJECTID | Id | Shape_Leng | Shape_Area | geometry | |
|---|---|---|---|---|---|
| 0 | 1 | 0 | 2000.000000 | 250000.000000 | POLYGON ((18162.733 31755.775, 17662.733 31755... |
| 1 | 2 | 0 | 2000.000000 | 250000.000000 | POLYGON ((23162.733 42755.775, 23662.733 42755... |
| 2 | 3 | 0 | 2000.000000 | 250000.000000 | POLYGON ((37662.733 40255.775, 37162.733 40255... |
| 3 | 4 | 0 | 2961.516312 | 226945.081982 | MULTIPOLYGON (((28662.733 48255.775, 28235.903... |
| 4 | 5 | 0 | 1974.913982 | 248059.096956 | POLYGON ((29162.733 28255.775, 29162.733 27755... |
gdf_grid = gpd.sjoin(gdf,grid_layer,op="within")
gdf_grid.head()
| CASE_ID | CLUSTER_ID | Onset_EYear | Onset_EWeek | Clean_Onset | ONSET_DATE | NEA_ONSET_ | NOTIFY_DAT | POSTAL | X output | ... | cell_id | friendly_n | Shape_Leng_left | Shape_Area_left | geometry | index_right | OBJECTID | Id | Shape_Leng_right | Shape_Area_right | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | DENF-2013-000025 | 0 | 2012 | 51 | 17/12/2012 | 17/12/2012 | 17/12/2012 | 31/12/2012 | 650536 | 18857.18765 | ... | 1334 | HONG KAH NORTH | 0.016238 | 0.000019 | POINT (18857.188 37512.049) | 2643 | 2644 | 0 | 2000.0 | 250000.0 |
| 95 | DENF-2013-000096 | 20130036 | 2013 | 1 | 5/1/2013 | 5/1/2013 | 5/1/2013 | 7/1/2013 | 650533 | 18705.23952 | ... | 1334 | HONG KAH NORTH | 0.016238 | 0.000019 | POINT (18705.240 37551.897) | 2643 | 2644 | 0 | 2000.0 | 250000.0 |
| 511 | DENF-2013-000530 | 20130036 | 2013 | 3 | 14/1/2013 | 16/1/2013 | 14/1/2013 | 19/1/2013 | 650532 | 18695.52780 | ... | 1334 | HONG KAH NORTH | 0.016238 | 0.000019 | POINT (18695.528 37593.481) | 2643 | 2644 | 0 | 2000.0 | 250000.0 |
| 773 | DENF-2013-001011 | 20130087 | 2013 | 4 | 26/1/2013 | 26/1/2013 | 26/1/2013 | 1/2/2013 | 650535 | 18772.05585 | ... | 1334 | HONG KAH NORTH | 0.016238 | 0.000019 | POINT (18772.056 37451.456) | 2643 | 2644 | 0 | 2000.0 | 250000.0 |
| 1018 | DENF-2013-001029 | 20130087 | 2013 | 5 | 27/1/2013 | 27/1/2013 | 27/1/2013 | 1/2/2013 | 650535 | 18772.05585 | ... | 1334 | HONG KAH NORTH | 0.016238 | 0.000019 | POINT (18772.056 37451.456) | 2643 | 2644 | 0 | 2000.0 | 250000.0 |
5 rows × 21 columns
grid_subset = grid_layer[grid_layer["OBJECTID"].isin(gdf_grid["OBJECTID"])]
grid_subset=grid_subset.to_crs(epsg ="4757")
grid_subset.to_crs(epsg = "3414")["geometry"].area/10**6
3 0.226945
4 0.248059
6 0.250000
8 0.250000
13 0.250000
...
3490 0.250000
3491 0.250000
3495 0.250000
3496 0.250000
3497 0.250000
Length: 1482, dtype: float64
import folium
from folium.plugins import HeatMap
map_folium = folium.Map([1.35255,103.82580],zoom_start=11.4,control_scale=True)
for _, r in grid_subset.iterrows():
#without simplifying the representation of each borough, the map might not be displayed
#sim_geo = gpd.GeoSeries(r['geometry'])
sim_geo = gpd.GeoSeries(r['geometry']).simplify(tolerance=0.001)
geo_j = sim_geo.to_json()
geo_j = folium.GeoJson(data=geo_j,style_function=lambda x: {'fillColor': 'orange'})
# folium.Popup(r['BoroName']).add_to(geo_j)
geo_j.add_to(map_folium)
map_folium